


//scatter plot with linear regression

*main first stage results gives the points for plotting the regression line

reghdfe ed_prescription propensity i.age_bin race_white female junior_enlisted college married longevity afqt_p , absorb(hym hd i.diag mos) cluster(provID)
predict predict_rate
local b = string(_b[propensity], "%12.3f")
local se = string(_se[propensity], "%12.3f")

*Residualize the Y variable

reghdfe ed_prescription i.age_bin race_white female junior_enlisted college married longevity afqt_p , absorb(hym hd i.diag mos) cluster(provID) resid(resid)


*Add back the mean for the residualized Y variable for ease of interpretation

replace resid = resid + _b[_cons]
collapse (mean) propensity resid predict_rate, by(provID year)

summ resid, d

*create figure

scatter resid propensity, mcolor(black) msize(medlarge) m(circle) ///
|| lfit predict_rate propensity, lcolor(black) lc(red) lw(.5) ///
title("", size(medlarge) color(black)) ///
ytitle("ED Prescriptions", size(medlarge)) ///
xtitle("Provider Propensity", size(medlarge)) ///
legend(off) graphregion(color(white)) 

graph export "$plotdir/resid_plot.png", as(png) replace





//local linear regression with confidence intervals

use "$data_out/provider_propensity.dta", clear


*residualize fe
reghdfe ed_prescription i.age_bin race_white female junior_enlisted college married longevity afqt_p , absorb(hym hd i.diag mos) cluster(provID) resid(resid)
replace resid = resid + _b[_cons]

lpoly resid propensity, nograph degree(1) gen(fs_x fs_y) n(100) se(se)

*store data
keep fs_x fs_y se 
drop if fs_x ==.
sort fs_x

save "$data_out/lpoly", replace


***load full dataset and create figure***
use provID year visitID propensity ed_prescription using "$data_out/provider_propensity.dta", clear
bys provID: egen raw_rate = mean(ed_prescription)
bys provID year (visitID): keep if _n ==1
append using "$data_out/lpoly"
gen upper = fs_y + 1.96*se
gen lower = fs_y - 1.96*se

summ propensity,d

local rmin `r(p1)'
local rmax  `r(p99)' 

display `rmin'
display `rmax'

tw histogram propensity if inrange(propensity, `rmin', `rmax'), color(gs10%30) percent yaxis(1) ///
|| line fs_y fs_x if inrange(fs_x, `rmin', `rmax'), lc(black) lw(0.6) yaxis(2) ///
|| line upper fs_x if inrange(fs_x, `rmin', `rmax'), lc(gs8) lw(0.3) yaxis(2) lp(dash) ///
|| line lower fs_x if inrange(fs_x, `rmin', `rmax'), lc(gs8) lw(0.3) yaxis(2) lp(dash) ///
title("") ///
ytitle("Fraction of Sample", size(medlarge) axis(1)) ///
ytitle("Prescription in ED", size(medlarge) axis(2)) yscale(range(0 .3) axis(2)) ///
xtitle("Provider Prescribing Intensity Instrument") ///
graphregion(color(white)) legend(off) ///
saving("local_linear.gph" , replace)

*graph save Graph "$plotdir/local_linear.gph", replace

graph export "$plotdir/local_linear.png", as(png) replace

//second figure with raw distribution\

summ propensity,d

local rmin `r(p1)'
local rmax  `r(p99)' 

summ raw_rate,d

local rmin2 `r(p1)'
local rmax2  `r(p99)' 

display `rmin'
display `rmax'

tw histogram raw_rate  , color(gs10%30) percent yaxis(1) ///
title("") ///
ytitle("Fraction of Sample", size(medlarge) axis(1)) ///
xtitle("Provider Prescribing Rate") ///
graphregion(color(white)) legend(off) ///
saving("local_linear_raw.gph" , replace)

*graph save Graph "$plotdir/local_linear_raw.gph", replace

graph export "$plotdir/local_linear_raw.png", as(png) replace




graph combine "local_linear.gph" "local_linear_raw.gph" , graphregion(color (white))  ycommon cols(2) ysize(2) 
   graph export "$plotdir/plot.png", replace


